Dynamical dimer method for the determination of transition states with ab initio 

molecular dynamics 
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A dynamical formulation of the dimer method for the determination of transition states is pre- 
sented. The method is suited for ab-initio molecular dynamics using the fictitious Lagrangian for- 
mulation. The method has been applied to the con-rotatory ring opening of chloro-cyclo-butadiene, 
an example, where the application of the drag method is problematic. 



I. INTRODUCTION 

The concept of the transition state has a fundamen- 
tal role for the prediction of rate constants of materi- 
als processes. The transition state is the lowest point 
on the energy barrier separating two metastable states 
representing the initial and final state of a process. The 
energy difference between the initial state and the transi- 
tion state is the activation energy. The curvatures of the 
potential energy surface at the initial state and the tran- 
sition state provides an estimate of entropic effects [1 [9]. 
A molecular dynamics simulation starting from the tran- 
sition state allows one to construct the dynamical tra- 
jectory of a reaction in the low-temperature limit. The 
extension of this latter approach to finite temperatures 
is straightforward [3|. 

The determination of transition states is however sub- 
stantially more difficult than the determination of min- 
ima of the total energy surface. A large number of meth- 
ods have been invented and refined over the last decades 
(see e.g. the Refs. fl H E9, M M 0, M M 

and the review articles [la. fl7j]). 

Many of these methods require the determination of 
the second derivatives of the total energy surface. While 
forces, that is first derivatives of the potential energy sur- 
face, are readily accessible from density functional calcu- 
lations, second derivatives are not, at least not without 
substantial effort. 

An overview of methods which allow to determine tran- 
sition states solely from force information can be found 
in [16f . The most simple approach is the drag method, 
where a one-dimensional constraint forces the system 
across a reaction barrier. Widely used are also "chain 
of states" algorithms such as the nudged elastic band 
(NEB) [19, [SOI or the strm S method [II El- 

Intermediate between these extremes is the dimer 
method. [U|. It couples two instances of the system, 
which we will call monomers. The monomers have a 
specified distance in configuration space. If the system 
follows the forces, after inverting of the force parallel to 
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the dimer, it evolves to a saddle point of first order, that 
is the transition state. 

The original dimer method [HI applied a two step 
algorithm. In a first step, the potential energy of the 
dimer-construct is minimized, excluding parallel motion, 
following a conjugate gradient approach. Then a large 
step is taken in the direction parallel to the dimer axis. 
Recently, an improvement has been suggested which min- 
imizes the number of gradient evaluations [19] • 

In the context of fictitious-Lagrangian a ppr oach to 
first-principles molecular dynamics (FPMD) [20], such a 
two-step procedure is not feasible, because the electronic 
wave functions must be able to follow the atomic mo- 
tion adiabatically. The fictitious Lagrangian approach 
requires a larger number of time steps. However the 
computational cost per time step are minimal. While 
a first-principles string approach has been published re- 
cently [211 ) . the dimer method has not been examined in 
the framework of ab initio molecular dynamics. 

In the present study, we introduce a dynamical formu- 
lation of the dimer method. Starting from an extended 
Lagrangian for the dimer in the configuration space with 
doubled dimensionality, we derive equations of motion, 
that, in the presence of dissipation, converge at transi- 
tion states of first order. An analysis of the trajectories 
of the dimer close to the stable points provides informa- 
tion on the stability and gives guidance for improving the 
performance of the method. The method has been im- 
plemented in the CP-PAW code and has been applied to 
the con-rotatory ring-opening of chloro-cyclo-butene, an 
example where the conventional drag method fails unless 
special precautions are taken. 

The paper is organized as follows: In section [IT] we de- 
fine the extended Lagrangian for the dimer motion and 
derive the equations of motion. In section IIIII we ana- 
lyze stability and the dynamics in the proximity of the 
stable points. Section Hvl is devoted to the discretization 
the equations of motion and the aspects of the actual 
implementation. In the final section [Vj the method is 
applied to the con-rotatory ring opening of cyclo-butene 
as practical example and compared to the drag method. 
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II. LAGRANGIAN AND EQUATIONS OF 
MOTION 



In this section we define the Lagrangian underlying 
the present work and derive the equations of motion. We 
begin by defining our notation and we introduce a con- 
venient set of coordinates: 

The dimer consists of two points in configurational 
space separated by a fixed distance. We use mass- 
weighted coordinates defined as 

x» := m'R,, (1) 

where the vectors Ri with i = 1,2 are the positions of 
the two monomers forming the dimer in configurational 
space. For a system with N atoms each vector has 3N 
dimensions. In addition, each monomer has its own set of 
electronic wave functions. The mass matrix m is diagonal 
and the masses of the respective atoms are located on the 
main diagonal. 

It is convenient to introduce a variable transformation 
into center-of-gravity coordinates Yi and relative coor- 
dinates Y 2 . 



Yr 



Xl +x 2 



Y 2 = xi - x 2 



(2) 
(3) 



Loosely speaking Yi describes the mean structure of the 
molecule, while Y 2 describes the orientation of the dimer. 

A dimer has three basic types of motion: motion par- 
allel to the dimer axis (vy), motion perpendicular to the 
axis (vj_) and rotation about the center of gravity (v D ). 
In the following, we divide the velocities into the con- 
tributions from these basic types of motion. The corre- 
sponding velocities are 



Y 2 = xi - x 2 
Y 2 ®Y 2 



d 2 

(xi - x 2 

Y 2 ®Y 2 



Yi 

(xi - x 2 ) (xi 



x 2 



1 



(f- 



2d 2 
Yi = 



Xl + x 2 



(4) 

(5) 
(G) 



The dot denotes the time derivative, and d is the so called 
'dimer distance' d = \/Y% = J (xi — x 2 ) 2 . The sym- 
bol ® refers to the outer product defined by (a ® b) c = 
a(b • c). Note, that Y2 Js Y2 is an operator that projects 
a vector onto the dimer axis. 

Now we define our theory by setting up a Lagrangian 
function in the relative and center-of-gravity coordinates 

£(Y 1 ,Y 2 ,Y 1 ,Y 2 ) = M o ^v 2 +M L 



M||V 



-V 



-V 



Yi 
Yi 



tY, 



-A[Y 



(7) 



The potential energy is described by the potential en- 
ergy surface V(R). The velocities v , vj_ and V| arc 
to be treated as functions of Yi and Y 2 and their time 
derivatives Yi and Y 2 as defined in Eqs. 1415(51 

We introduced three scale factors M , My and Ml, 
that allow us to scale the masses for the three types of 
motions independently, that will later be used to accel- 
erate convergence. More importantly, by choosing a neg- 
ative value of M|| , the motion along the dimer direction 
is inverted, so that the dimer climbs up the barrier. 

If the scale factors are equal to one, the conventional 
kinetic energy is recovered, that is 



4 V ° 



vf + vi = ^RimRi + iR 2 mR 2 (8) 



Using the method of Lagrange multipliers, a term 
has been introduced to describe the dimer-distance con- 
straint 



(9) 



which ensures that the dimer distance in mass-weighted 
coordinates is equal to d. The corresponding Lagrange 
multiplier is denoted by A. 

In addition to the Lagrangian we also define a 
Rayleigh's dissipation function T> 



V 



1 



loM -wi + 7||M||Vjj + 7±M ± vi . (10) 



which allows us to introduce dissipation in a consistent 
manner. 

The Euler-Lagrange equation for the Lagrangian of 
Eq.[7]and the Rayleigh's Dissipation function from Eq.fTUl 
are obtained from 



d dC 



dC 



dV 



dt 3Y, n dY i>n dY itn 



0, 



(11) 



where the index i € {1, 2} labels the two monomers, and 
n 6 {l,3iV} labels the coordinate in configuration space. 
The resulting equations of motion have the form 



" Y 2 ®Y 2 
d 2 



My (Y 1 + 7 ||Y 1 
1 



-m"5 (Fi +F 2 ) 



(M { \ - Mi 



d 
dt 



Y 2 



d 2 



Y x = (12) 



and 



M Y 2 -m-5(F 1 -F 2 
Y 2 Yi 



+4 (M ± - Mi 



d 2 



4AY 2 

Yi + 7oM Y 2 = (13) 
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Here we have used the forces acting on the monomers 
defined as 



F, := - V 



R, 



V 



(14) 



In equation (fT3"|) we absorbed all terms which lead to a 
force parallel to the dimer axis into the constraint force 
by redefining the Lagrange multiplier. Thus the variable 
A used in Eq. Q2] differs from the Lagrange multiplier A 
in Eq. 

Furthermore, because the dimer system moves on the 
hyperplane with constant dimer distance, we simplified 
the final equations by using 



Yn = d 2 



and 



—~Y Z 

at 2 



o. 



(15) 



(16) 



Equation (fT2"|) describes the motion of the center of grav- 
ity of the dimer, which is related to the structure of the 
molecule. Equation (|13|) describes the orientational mo- 
tion of the dimer. 



III. LOCAL STABILITY ANALYSIS 

For M = Mm = Ml = 1, the Lagrange function ([7]) 
describes a physical system of two masspoints moving in 
a potential V under the influence of the constraint force 
that keeps the dimer distance invariant. With positive 
friction factors 7 , 7|i and j± the dimer will come to rest 
with the center-of-gravity coordinate next to a local min- 
imum. The dimer axis will be aligned nearly parallel to 
the lowest vibrational eigenmode. The above statements 
are exactly fulfilled in the limit of vanishing dimer dis- 
tance. 

If we choose a negative value of Mm , the motion will be- 
come unstable near local minima. Instead, the dimer will 
be attracted by transition states of first order. A transi- 
tion state of first order is characterized by the presence 
of exactly one eigenmode with an imaginary frequency. 

These properties of the dynamics have been derived 
from the following local stability analysis. First we de- 
termine the stationary points for the dimer dynamics. 
Then we investigate the dynamics in the neighborhood 
of those stationary points. 

For this analysis we can replace the potential energy 
surface by its truncated Taylor expansion at a given point 
R(°) — m"5x' ', namely 



V(R) = W 0) -F(°>R 



R(°))m5Dm5(R-R(°)) (17) 



with the dynamical matrix D defined as 



D 



1 



d 2 V 



1,3 



1 



R(0) V 7 



(18) 



From Eq. [17] and Eq. [T] we obtain the forces 

F(x) =P (0) -msD(x-xf°)) (19) 

A. Stationary points 

If we insert the condition for a stationary point Y^ = 
Y s ; = in the equation of motion Eq. [T2l we find that 
the forces for the two monomers at the stationary point 
are antiparallel and of the same magnitude. Using the 
Taylor expansion Eq.[I]J]m the equation of motion Eq. ll2[ 
we find that the center-of-gravity of the dimer Yi lies 
at a stationary point of the potential, if the dimer is 
stationary, that is 

F(Yi) = (20) 
Similarly we obtain from Eq. Q1J] and Eq. [TH] 



DY 2 = -4AY 2 



(21) 



Thus the dimer axis points along an eigenvector of the 
dynamical matrix. 

In conclusion we find that the dimer dynamics is sta- 
tionary, when its center of gravity lies at an extremum or 
a saddle point of the potential energy surface and when, 
in addition, the dimer axis points along one of the vibra- 
tional eigenmodes. 



B. Linearized equation of motion 

Without loss of generality, we choose in the following 
a coordinate system for which the stationary point of the 
potential energy surface lies at the origin. Secondly, we 
denote the eigenvectors of the dynamical matrix as e^, 
with e,;ej = Sij, and the corresponding eigenvalues with 
to 2 . The eigenvector, which is parallel to the dimer axis is 
denoted by e and the corresponding eigenvalue is denoted 
by Co 2 . Thus, with Eq. [2U we can identify the Lagrange 
multiplier as A = — ju> 2 . 

Linearization of the equations of motion Eqs. [12] and 
[H about Y 1 = and Y 2 = ed with = yields the 
following equations for the deviation SYi(t) and SY2(t) 
from the stationary point 



(1 - e <8> e) Mi ( <5Yi + -/^SY 1 

+ (e®e)M|| (SYi +7||<5Yi) +T>SY 1 = (22) 
MJY 2 + BSY 2 + Cj 2 5Y 2 + M olo 8Y 2 = (23) 

We project the second equation ,Eq.[221 onto the eigen- 
vectors ei and obtain: 

M (e 4 <5Y 2 ) + (lj 2 - w 2 )(e j( 5Y 2 ) + M y (eiSY 2 ) = 0(24) 

Note, that the motion along e is simple due to the dis- 
tance constraint. The frequency of the variable e^Y 2 is 



We see that, for a positive mass M , the dynamics is 
stable only, if w 2 is the lowest eigenvalue of the dynamical 
matrix. Thus the dimer axis will always orient along the 
eigenvector with the lowest eigenvalue. 

Now, we project Eq. [22] onto the eigenvector e of the 
dynamical matrix 

M||(e<5Yi) + cj 2 (e(5Yi) + M||7||(e(SYi) = (26) 

which yields the translational motion of the dimer along 
the dimer axis. This equation will be responsible for 
the ascent to a saddle point. The eigenfrequency of the 
variable eYi is 



wii = ±- 



(27) 



For a negative mass Mm , which is the choice for a tran- 
sition state search, Eq. [55] is stable only if the dimer is 
oriented along an unstable vibrational mode of the po- 
tential energy surface, that is with the dimer oriented 
along an eigenvector e with an imaginary frequency. 

Now, we project Eq. [5H onto the eigenvectors with 
e,- _L e and we obtain 



MxMYr) 



2 (e l SY 1 ) + M^ ± (e l SY 1 ) = (28) 



which represents the translational motion perpendicular 
to the dimer axis. The translational motion is related to 
an optimization of the atomic structure. The frequency 
of the variable e^Yi with _L e is 



FIG. 1: Sketch of the eigenvalue spectrum of the dynamical 
matrix at the saddle point of first order with the definition of 

Ldmin and UJmax- 



Thus, the stability limit requires us to choose the 
masses as 



Q 2 A 2 



Mi, < 



. .2 A 2 
^ max 
± > ^ 



(2 _ 2 U2 
^ j ^ v^max ^min) 



(30) 
(31) 
(32) 



(29) 



In order to minimize the number of time steps, it is de- 
sirable to choose the masses close to these limits. 

Let us now obtain an a priori estimate of these lim- 
its: (1) For practical purposes, we can make the assump- 
tion that Lu m i n << oJmax an d thus we set 0J m in to zero. 
(2) The highest vibrational mode is that of H 2 , with a 
frequency of 4400 cm -1 . Frequencies for bond-stretch 
vibrations are of the order 1000 cm -1 . These numbers 
are reasonable upper estimates for uj max . (3) The abso- 
lute value of the imaginary frequency at the saddle point 
is typically in the range of the high-lying real-frequency 
modes, that is comparable to ui m ax- Hence uj 2 ~ —^ max - 
With these assumptions we obtain the following, recom- 
mended values for the masses 



From Eqs. [25] [27] and [29] we find that the dimer with 
positive M and M± and negative Mm will only come to 
rest at transition states of first order, that is, if Co 2 < 
and if all other eigenvalues are positive. In that case the 
dimer axis points along the eigenvector corresponding to 
the imaginary frequency, that is, across the saddle point. 



C. Masses and Frictions 

We will later solve the equations of motion [T2l [TBI with 
the Verlet algorithm. The Verlet algorithm for a har- 
monic oscillator with a frequency u> becomes unstable, if 
the time step A used for the discretization of the equa- 
tion of motion is smaller than 2/u. The frequency of the 
discrctized motion is accurate to within 1 % if wA < -| . 
Thus, it is important to understand the vibrational spec- 
trum, shown schematically in Fig.[TJ of the motion in the 
potential. 



My 

Ml 
M 



— kA 2 

kA 2 

2kA 2 



(33) 
(34) 
(35) 



where n = (6-10 3 a.u.) 2 is a recommended constant 
based on the H2 vibration. 

Let us now turn our attention to the friction values. 
The goal is to reach an optimum convergence at the sta- 
tionary point of the dimer. For a damped harmonic os- 
cillator 



mi 



-muf^x 



rwyx 



(36) 



the fastest decay rate as function of the friction is ob- 
tained for critical damping, that is in between the oscil- 
latory and the over-damped regime. Critical damping is 
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obtained for a friction of 7 = 2ujq. Thus we choose 

< 7-1 < 2 




9 / ^rnin 

V Mi 



Mi 



w 2 . -d» 2 

mm 



Mo 



< 7o < 2i 



tut 



Mo 



(37) 
(38) 
(39) 



Note, that a too high friction freezes out those modes in 
the over-damped regime. Thus it is usually better to use 
a friction near the lower bound of the sensible regime. 



IV. NUMERICAL INTEGRATION OF THE 
EQUATIONS OF MOTION 

A. Discretization of the equations of motion 



The equations of motion Eqs. 1121131 are nonlinear in 
the velocities. This leads to a nonlinear equation for the 
discrctized equation of motion. We tackle the problem 
by iterating on the nonlinear terms in the velocities. 



We set up the discretized equations of motion, while 
treating the terms nonlinear in the velocities as an ab- 
stract force. 



G 2 



(My - Mi) 



d_ (Y 2 ®Y 2 
dt 



d 2 



Yi (40) 
(41) 



The equations of motion are discretized according to 
the Verlet algorithm: 



y 
v 



v{+) - y{-) 

2A 

y(+)-2y(0)+y(-) 



(42) 
(43) 



where A is the discretization time step, and where we 
used the short-hand notation 



y{+) = y{t + &) 

dy(0) = y(t) 
y(-) = y(t-A) 



We obtain from Eq. [T^] 



(44) 
(45) 
(46) 



1 - 



+ 



Y 2 ® Y 2 \ M ± 



d 2 



Y 2 \ M, 



d 2 



^{(l + ox)Yx(+)- 
^{(l + «ll)Yi(+) 



2Y 1 (0) + (l-o L )Y 1 (-)} 
^YxtOj+fl-aiiJYiH} 



1 



m 2 (p 1 + F 2 ) - Gi 



(47) 



where aj_ = and an — . In the projectors we The result is 
have dropped the argument of Y 2 (0). 

Eq. d7] can be resolved for Yi(+) by multiplication 
with 



1 



Y 2 ®Y 2 \ A 



d 2 J M ± (l + a ± ) 



v2 



. Y 2 ® Y 2 \ A 2 
+ 1— ^— J M ||(l + a||) (48) 
I 



Yl(+) _ ^.^{^(oj.^h^^,^^).^} 



Y 2 ®Y 2 
d 2 



(t^-Y^O) - i-^Y x (-) + (im-*(F! + P 2 ) + G 1 ) f } (49) 
l-l + 0|| l + a || V 2 /M||(l + a||)J 

r 



Similarly we discretize the equation of motion for the dimer distance, Eq. [T5] Following Ryckaert et al [22], we 
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first propagate without the force of constraint to obtain 

a n 



Y 2 = -^-Y 2 (0)-i- 

1 + a H 

fm-i(Fi-P 2 ) 



-Y a (-) 



A 2 



, , , \~ (5°) 
M (l + a ) 

with a — The new vector Y 2 (+) is related to Y 2 

by the constraint force and the Lagrange parameter A 
according to 

Y 2 (+)=Y 2 -4AY 2 (0) f (51) 

Mo(l + do) 

The Lagrange parameter is then adjusted so that the 
constraint Y 2 2 (+) = d 2 , namely a given dimer length, is 
satisfied. 

In the first iteration we estimate this forces Gi and 
G 2 from the previous iterations or we set them to zero. 
Then we propagate the positions, which provides us with 
a better estimate for the nonlinear terms. This loop is 
then iterated to convergence. 



B. Restricting the orientation of the dimer 

It will be important to limit the degrees of freedom, 
in which the two monomers may differ. This implies re- 
stricting the orientation of the dimer to a space with 
lower dimensionality. The reason is to avoid that the 
dimer converges at irrelevant saddle-points, which are not 
of interest. This is a common problem for complex sys- 
tems, that exhibit many minima and saddle points. This 
restriction is accomplished easily, by setting the corre- 
sponding components of Y 2 in Eq. [50] to zero. 



V. APPLICATION: CONROTATORY RING 
OPENING OF l-CHLORO-2-CYCLOBUTENE 

In order to demonstrate the performance of the dimer 
method described above, we have chosen a prototypical 
system for tests of transition state searches, namely the 
ring opening of cyclobutene. The isomerization of cy- 
clobutene to cis-butadiene is the prototypical example of 
concerted stereospecific reactions. The underlying pro- 
cesses have been studied and discussed extensively by 
several groups [H, [24], [25|, [26[ (and references therein). 
The potential energy surfaces of the con- as well as of 
the disrotatory isomerization mechanisms exhibit prin- 
cipal structures, which render the determination of the 
corresponding transition states complicated [24j. 

In order to avoid special effects due to the high sym- 
metry of cyclobutene, that would be untypical for other 
molecules, we explored l-chloro-2-cyclobutene. This 
molecule and the two relevant reaction products are 
shown in Figure [2j 

Computational details of our calculations are given in 
appendix [X] The coordinates of the structures are sup- 
plied with the supplementary material. 
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FIG. 2: Initial state, l-chloro-2-cyclobutene (A), and final 
state, l-chloro-buta-l,3-diene (B), for the conrotatory ring 
opening of chlorocyclobutene. Also shown is the product 
of the disrotatory ring opening, trans- l-chloro-buta-l,3-diene 
(C). 



A. Drag Calculations 

In this section we investigate the reaction using the 
drag method. Because of its simplicity, the drag method 
is widely used for transition state search. However, the 
drag method fails for a number of systems. In such cases 
the energy exhibits a hysteresis effect, that is, different 
energy profiles are obtained when dragging the reaction 
coordinate in the opposite directions. The ring opening 
of cyclobutene is one example where the drag method 
exhibits a hysteresis effect. 

In the drag method one specifies a one-dimensional 
reaction coordinate. Then one maps the energy profile 
along the reaction coordinate and determines the highest 
point as the transition state. Each point along the en- 
ergy profile corresponds to the local energy minimum on 
a hyperplane perpendicular to the reaction coordinate. 
In our case we selected the reaction coordinate by the 
difference between initial and final state. A hyperplane 
with a specified value c of the reaction coordinate is given 
by 



(Kb - Ra) (R - Ra) 
(Rb — Ra) 2 



(52) 



Hence the value of the reaction coordinate is zero for the 
initial state (A) and one for the final state (B). 

In addition to the reaction coordinate we imposed six 
additional constraints to avoid translations and rotations 
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of the molecule. We have chosen 

i(R(C 2 ) + R(C 3 )) = (53) 

z(d) = z(C 4 ) = I {x{C x ) + x(C 4 )) = (54) 

Where the coordinates correspond to the atoms given in 
parenthesis. The notation follows Figure [H 

When the energy profile is determined by varying the 
reaction coordinate in small steps, once from zero to one 
and then from one to zero, we obtain the hysteresis shown 
in Fig. [3l The highest point of the energy profile is not 
the transition state. Instead, the highest point of each 
branch is an upper bound for the activation energy, while 
the crossing of the two branches is a lower bound. 



1 • 1 ■ 

: / 


i ■ 1 ■ 1 ■ r 







0.0 0.2 0.4 0.6 0.8 1.0 



c 



FIG. 3: Potential energy relative to the energy of final state 
(B) over reaction coordinate from a forward and backward 
calculation using the drag method. The true transition state 
(determined using the dynamical dimer method) is denoted 
by a circle. 

In order to show the underlying reason for the hys- 
teresis, we show in Fig. 2] the total energy surface in a 
two-dimensional hypersurface. The two axes are the re- 
action coordinate and the coordinate defined by the two 
isoenergetic structures at the crossing of the two branches 
of the energy profile in Fig. [3J For each point in the two- 
dimensional plot the energy is at a local minimum of the 
(3./V — 2)-dimensional hypersurface. 

We see that the drag method in the forward direction 
leads into a side valley that leads further onto a ridge. 
The point where this path becomes unstable, and from 
where it leads down into the valley of the product state, 
lies past the transition state. In the backward direction 
a similar instability occurs, even though closer to the 
reaction coordinate of the transition state. 

Note that the two configuration A and , that mark 
the crossing of the two branches in the energy profile 
shown in Fig. [3J lie far from the transition state. For 
further information see Table |TJ 




0.2 0.4 0.6 0.8 1.0 



c 

FIG. 4: Sketch of the potential energy surface of the con- 
rotatory isomerization of l-chloro-2-cyclo-butene to 1-chloro- 
buta-l,3-diene. The dashed lines trace the paths obtained 
from drag calculations for the transitions A — > B and B — > A, 
respectively. The crosses denote the structures A* and . 
TS corresponds to the transition state. 



B. Dynamical dimer calculations 

In this section we describe technical issues for dy- 
namical dimer calculations and demonstrate their per- 
formance. 

We use the fictitious Lagrangian formulation of ab- 
initio molecular dynamics. This implies that wave func- 
tion coefficients and nuclei obey Newton's equation of 
motion, to which we added a friction, which allows to 
quench the system. 

Special attention should be given to the wave function 
dynamics. The wave- function cloud tied to the atoms re- 
sults in an increased effective mass of the nuclei. Further- 
more, a friction applied to the wave function dynamics 
acts like an effective friction on the nuclear motion. For 
a normal ground-state search, this is not problematic. 
In our case however, the mass for the motion parallel to 
the dimer direction is inverted. Thus the dimer acceler- 
ates opposite to the direction of the force. Thus a friction 
acting on the dimer via the electrons leads to an velocity- 
dependent acceleration of the dimer, that may cause an 
instability of the dimer motion. 

The detrimental effect of the wave function motion can 
be avoided by choosing a sufficiently small mass for the 
wave function dynamics or by artificially increasing the 
atomic masses. In our study, we used a mass of 50 u for 
all atoms. The wave functions have been kept close to 
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the ground state by applying a friction that dissipates 
2 % of the kinetic energy in each time step. 

The dimension-less masses for the dimer motion have 
been chosen to 

My = -1.00 (55) 
Mi = +1.00 (56) 
M = +0.25 . (57) 

The small value for M Q has been chosen to speed up the 
reorientation of the dimer. 

In order to avoid instabilities, we found it useful to 
introduce an upper limit to the kinetic energy individu- 
ally for the rotational, perpendicular and parallel motion. 
These limits were enforced by adjusting the frictions ac- 
cordingly. The limit for the kinetic energy Ekin,max is 
translated into an upper "temperature" T max according 
to 

~^Qkj^T rnax — Eki nrnax (^8) 

where g is the number of degrees of freedom in the cor- 
responding motion type, kg is Boltzmann's constant. 

Within the limits of the enforced maximum kinetic en- 
ergies, we adapted the friction dynamically in each step 
to come close to the limit of critical damping. The latter 
results in the best possible convergence behavior. 

Furthermore, the friction is increased to a higher value 
specified by the user, if the dimer moves away from the 
stable point. The decision to increase the friction has 
been made on the basis of the actual forces and velocities. 
In our simulation we have chosen this friction so that the 
energy dissipated per time step corresponds to 20 % of 
the kinetic energy of the corresponding type of motion. 

When starting the dimer simulation, it is important 
to optimize the dimer orientation, i.e. Yi, first, before 
the mean dimer configuration Y2 changes appreciably. 
This is because the mean dimer configuration only ap- 
proaches the transition state, if the dimer orientation is 
sufficiently aligned along the unstable vibrational mode 
of the transition state. 

The dimer orientation can be optimized in different 
ways. (1) Starting from initial dimer coordinates, the 
mean dimer configuration, Yi, is constrained to the ini- 
tial value, while the orientation, Y2, is fully optimized. 
Only after this optimization also the mean configuration 
is allowed to move. (2) Starting from one monomer con- 
figuration, the second monomer is constructed from the 
first by letting it follow the forces until the desired dimer 
length is obtained. We will refer to this technique as 
the "growing-dimer technique" (3) Starting from a dimer 
with the monomers being identical to the two metastable 
minima, the dimer length is successively shortened, while 
the dimer position is optimized. Results for all three op- 
timization strategies can be found in Table HI 

We found that the third strategy suffers from the fact 
that additional effort is required to contract the dimer 
length to the desired value. It also suffers from an insta- 
bility caused by strong anharmonic effects. 



The growing-dimer technique has the advantage that 
the second monomer is constructed dynamically from the 
first, without the need of an independent optimization of 
the wave functions. In the growing dimer technique the 
dimer orientation is optimized automatically as soon as 
the targeted dimer length is reached. In the first strategy 
the orientation is obtained during the first iterations with 
a constrained mean dimer configuration. 

We will discuss here the first strategy. The initial 
structure of the dimer is given by a mean dimer config- 
uration midway between the two metastable states (A) 
and (B). The dimer length is chosen such that the dis- 
tance of the two monomers in coordinate space is 0.33 A. 

During the first 700 time steps, the dimer orientation 
has been optimized. Here the mean dimer configuration 
has not been constrained, but the maximum temperature 
for the perpendicular motion has been kept at a small 
value of 10 K. For the parallel and rotational motion, 
the maximum temperature has been set to 500 K. After 
the first 700 time steps the maximum temperature for 
the perpendicular motion has been increased to 500 K. 
In Figs. 15161 we show the convergence of the forces for 
the parallel, rotational and perpendicular motion. Good 
convergence is reached after 2000 time steps. The esti- 
mates for the transition state energy and the deviation of 
the mean dimer configuration from the transition state 
is given in Table fl] 




'o 1000 2000 3000 

Simulation Step 



FIG. 5: Forces over time-step for a calculation following op- 
timzation strategy (1) described in the text. The dashed and 
dash-dotted line show the forces parallel to the dimer axis for 
image one and two, respectively. The full line corresponds 
to the resulting force acting on the center of gravity of the 
dimer. The logarithmic scale of the inset provides detailed 
information for the latter. 



VI. SUMMARY 

A formulation of the dimer method for searching tran- 
sition states is presented which can be used in ab-initio 
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FIG. 6: Rotational (full line) and perpendicular (dashed line) 
part of the forces acting on the dimer over time-steps for a 
calculation following optimzation strategy (1) described in the 
text. 



TABLE I: Predicted transition-state energy (relative to final 
state (B)) and deviation from the exact transition state for 
the con-rotatory isomerization of l-chloro-2-cyclo-butene to 
l-chloro-buta-l,3-diene. The results include three different 
dynamical dimer optimization strategies ((1), (2) and (3)) 
described in the text and the forward as well as the backward 
drag calculation shown in Fig. [3] 

Predicted Energy (eV) |R TS - R| (A) 



Transition State" 


2.110 




(1) 






1000 steps 


2.110 


1.00 


1500 steps 


2.111 


0.85 


2000 steps 


2.111 


0.71 


(2) 






1000 steps 


2.258 


0.62 


1500 steps 


2.116 


0.37 


2000 steps 


2.110 


0.36 


(3) 






1000 steps 


2.133 


1.58 


1500 steps 


2.112 


0.85 


2000 steps 


2.110 


0.85 


TS A -> B 


3.461 


2.57 


TS B -» A 


2.120 


0.56 



APPENDIX A: COMPUTATIONAL DETAILS 

We performed density-functional calculations [13, 
based on the projector augmented wave (PAW) 
methodJH HI- The gradient-corrected PBE[U func- 
tional was used for exchange and correlation. The PAW 
method is a frozen-core all-electron method. Like other 
plane-wave based methods, the PAW method leads to 
the occurrence of artificial periodic images of the struc- 
tures. This effect was avoided by explicit subtraction 
of the electrostatic interaction between them. [32| Wave 
function overlap was avoided by choosing the unit cell 
large enough to keep a distance of more than 6 A be- 
tween atoms belonging to different periodic images. We 
used a plane wave cutoff of 30 Ry for the auxiliary wave 
functions of the PAW method. The following sets of pro- 
jector functions were employed, CI 2s2pld, C 2s2pld, H 
2slp, which provides the number of projector functions 
per angular momentum magnetic quantum number m in 
each main angular momentum channel I. 

Atomic structures were optimized by damped Car- 
Parrinello[33] molecular dynamics. We used a time-step 
of 10 a. u. (2.5 fs) for all except the dynamical dimer cal- 
culations. See section IV Bl for further details. The con- 
vergence was tested by monitoring if the total energy 
change remains below 10 -5 Hartree during a simulation 
of 500 time steps. During the simulation for the conver- 
gence test, no friction was applied to the atomic motion 
and the friction on the wave function dynamics was cho- 
sen sufficiently low to avoid a noticeable effect on the 
atomic motion. 



a Determined by the use of a well converged dynamical dimer cal- 
culation with final dimer distance of 0.125 A. 



molecular dynamics simulations using a fictitious La- 
grangian. The dimer method is successful in cases where 
the widely used drag method fails. We investigate the 
dynamics close to the stable points of the dimer analyti- 
cally. In addition we demonstrate the performance of our 
implementation using a practical example typically used 
as test case for transition state search algorithms, namely 
the con-rotatory ring opening of (chloro-) cyclo-butene. 
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